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A  BOUNDARY- FITTED  COORDINATE  CODE  FOR  GENERAL  TWO-DIMENSIONAL 
REGIONS  WITH  OBSTACLES  AND  BOUNDARY  INTRUSIONS 


INTRODUCTION 

The  use  of  numerically  generated  boundary-fitted  curvilinear  coor¬ 
dinate  systems  as  the  basis  for  numerical  solution  of  partial  differential 
equations  on  arbitrary  regions  is  now  well  established.  A  comprehensive 
survey  of  the  generation  and  use  of  these  coordinate  systems  has  recently 
appeared,  Ref.  [1],  and  the  proceedings  of  a  recent  symposium  devoted 
to  this  area.  Ref.  [2],  cover  the  basic  techniques  involved,  as  well  as 
applications  in  many  areas. 

Such  coordinate  systems  have  the  property  that  some  coordinate  line 
is  coincident  with  each  segment  of  the  boundary  in  the  physical  region, 
so  that  the  complication  of  boundary  shape  is  effectively  removed  from 
the  problem.  In  the  past  decade  the  numerical  generation  of  curvilinear 
coordinate  systems  has  provided  the  key  to  the  development  of  finite 
difference  solutions  of  partial  differential  equations  on  regions  with 
arbitrarily  shaped  boundaries.  Although  much  of  the  impetus  for  these 
developments  has  come  from  fluid  dynamics,  the  techniques  are  equally 
applicable  to  heat  transfer,  electromagnetics,  structures,  and  all  other 
areas  involving  field  solutions. 

With  coordinate  systems  generated  to  maintain  coordinate  lines 
(surfaces  in  3D)  coincident  with  the  boundaries,  finite  difference  codes 
can  be  written  which  are  applicable  to  general  configurations  without 
the  need  of  special  procedures  at  the  boundaries.  Even  when  the  bound¬ 
aries  are  in  motion,  the  use  of  such  coordinate  systems  allows  all 
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computation  to  be  done  on  a  fixed  grid  with  a  uniform  square  mesh  in  the 
transformed  plane.  This  greatly  simplifies  the  coding,  particularly 
with  regard  to  boundary  conditions,  which  can  now  be  represented  without 
need  of  interpolation.  It  is  also  possible  to  distribute  the  curvilinear 
coordinate  lines  in  the  physical  plane  with  concentration  of  lines  in  regions 
of  high  gradients  while  maintaining  the  square  grid  in  the  transformed 
(computational)  plane. 

With  such  systems,  the  grid  points  may  be  thought  of  as  a  finite 
set  of  observers  of  the  physical  solution,  stationed  so  as  to  be  most 
effective  in  covering  all  of  the  action  on  the  field.  The  structure  of 
an  intersecting  net  of  families  of  coordinate  lines  allows  the  observers 
to  be  readily  identified  in  relation  to  each  other.  This  results  in 
much  more  simple  coding  than  would  the  use  of  a  triangular  structure 
or  a  random  distribution  of  points.  The  grid  generation  system  provides 
some  influence  of  each  observer  on  the  others  so  that  when  one  moves 
to  get  into  a  better  position,  its  neighbors  will  follow  in  order  to 
maintain  smooth  coverage  of  the  field.  The  curvilinear  coordinate  system 
thus  should  cover  the  field,  with  coordinate  lines  (surfaces)  coincident 
with  all  boundaries.  The  distribution  of  lines  should  be  smooth,  with 
concentration  in  regions  of  high  gradient. 

Numerical  solutions  of  partial  differential  equations  are  done  on 
the  curvilinear  coordinate  system  by  first  transforming  all  partial 
derivatives  (or  Integrals)  analytically  so  that  the  curvilinear  coordinates, 
rather  than  the  physical  coordinates,  become  the  independent  variables. 

Normal  and  tangential  derivatives  at  boundaries  are  similarly  transformed. 
(These  transformation  relations  are  given  in  Ref.  [3].)  The  result  is  a 
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set  of  partial  differential  equations  and  boundary  conditions  in  which 
all  derivatives  (and  integrals)  are  with  respect  to  the  curvilinear  coor¬ 
dinates.  These  equations  may  then  be  expressed  as  difference  equations 
on  the  square  grid  that  is  inherent  in  the  transformed  plane.  There  is 
thus  no  need  for  Interpolation  regardless  of  the  shape  of  the  boundaries 
or  the  distribution  of  the  curvilinear  coordinate  lines  in  the  field. 

The  present  report  concerns  a  code  for  the  generation  of  boundary- 
fitted  coordinate  systems  for  general  2D  regions  with  boundaries  of  ar¬ 
bitrary  shape  and  with  internal  obstacles  and  boundary  intrusions,  arbi¬ 
trary  in  shape  and  number.  The  code,  referred  to  as  WESC0R,  is  described 
and  instructions  for  input  and  use  are  given.  Examples  of  the  applica¬ 
tion  of  this  code  are  given  in  Ref.  [4]- [6].  The  coordinate  system  is 
generated  from  the  numerical  solution  of  a  system  of  elliptic  partial 
differential  equations  with  provision  for  controlling  the  spacing  of  the 
coordinate  lines  in  the  field.  The  transformed  (computational)  region 
is  rectangular  with  the  obstacles  and  intrusions  transformed  to  slits 
and/or  slabs.  (This  type  of  transformed  configuration  and  its  use  are 
discussed  in  Ref.  [3].)  A  small  code  to  distribute  points  on  various 
fundamental  curves  with  exponential  concentration  is  also  described. 

This  front-end  code  can  be  used  to  construct  boundary  point  distributions 
for  input  to  the  coordinate  code.  A  plot  code  for  the  coordinate  system 
is  also  included.  The  boundary-fitted  coordinate  systems  generated  by 
this  code  may  be  used  as  a  basis  for  the  numerical  solution  of  partial 
differential  equations  for  any  physical  problem  of  interest. 

Tne  elliptic  generation  system  is  discussed  in  Part  A,  and  the  op¬ 
eration  and  use  of  the  codes  are  covered  in  Part  B. 


PART  A 


ELLIPTIC  GENERATION  SYSTEM 


ELLIPTIC  GENERATION  SYSTEM 

The  generation  of  boundary-fitted  coordinates  from  elliptic  systems 
and  the  use  thereof  in  the  numerical  solution  of  the  Navier-Stokes  e- 
quations  is  surveyed  in  Ref.  [1],  The  foundations  of  elliptic  generation 
systems  are  discussed  in  detail  in  Ref.  [7],  and  basic  configurations  of 
the  transformed  plane  are  covered  in  Ref.  [3).  The  discussion  in  this 
section  is  an  introduction  to  the  subject  given  by  Johnson  and  Thompson 
in  Ref.  [5]  and  is  incorporated  here  for  convenience. 

Basic  Ideas 

Suppose  one  is  interested  in  solving  a  differential  system  involving 
two  concentric  circles,  such  as  shown  in  Fig.  1,  where  r  =  constant  “ 
on  the  inner  circle  and  r  =  constant  =  on  the  outer  circle,  and  0 
varies  monotonically  over  the  same  range  over  both  the  inner  and  outer 
boundaries,  i.e.,  0°  to  360°. 

A  cylindrical  coordinate  system  is  the  obvious  choice  since  a  coor¬ 
dinate  line,  i.e.,  a  line  of  constant  radius,  coincides  with  each  boundary. 
If  one  now  pulls  the  interior  regions  between  the  two  circles 


Figure  1.  Transformation  of  domain  between 
concentric  cylinders 
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apart  at  8  ■  0°  (or  0  •  360°)  and  folds  outward,  it  is  easy  to  visualize 
the  region  becoming  the  rectangular  region  D£.  Likewise,  it  should 
be  obvious  that  the  right  and  left  sides  of  the  rectangle  are  reentrant 
boundaries  since  0*0°  and  0  ■  360°  are  coincident  in  region  D^.  If 
one  computes  a  derivative  in  the  cylindrical  system  at  0*  0°,  values 
at  the  points  marked  x  and  o  on  both  sides  might  be  used.  Thus,  these 
same  points,  as  shown  in  the  rectangular  region,  would  be  used  for  a 
similar  derivative  in  region  D2.  This  is  the  reason  for  calling  these 
boundaries  reentrant  boundaries.  As  shown,  the  boundary  of  the  inner 
circle  becomes  the  bottom  of  the  rectangular  region  while  the  boundary 
of  the  outer  circle  becomes  the  top. 

The  general  boundary-fitted  system  is  completely  analogous  to  the 
system  discussed  above.  In  Fig.  2  the  curvilinear  coordinate,  rj ,  is 
defined  to  be  constant  on  the  inner  boundary  in  the  same  way  that  the 
curvilinear  coordinate,  r,  is  defined  to  be  constant  on  the  inner  circle 
in  the  cylindrical  coordinate  system.  Similarly,  n  Is  defined  to  be 
constant  at  a  different  value  on  the  outer  boundary.  The  other  curvi¬ 
linear  coordinate,  5,  is  defined  to  vary  monotonically  over  the  same 
range  on  both  the  inner  and  outer  boundaries,  as  the  curvilinear  coordi¬ 
nate,  0,  varies  from  0  to  2ir  around  both  the  inner  and  outer  circles  in 
cylindrical  coordinates.  It  would  be  just  as  meaningless  to  have  a  dif¬ 
ferent  range  for  £  on  the  inner  and  outer  boundaries  as  it  would  be  to 
have  0  Increase  by  something  other  than  2tt  around  one  of  the  circles  in 
cylindrical  coordinates.  It  is  this  fact  that  £  has  the  same  range  on 
both  boundaries  that  causes  the  transformed  field  to  be  rectangular. 

Note  that  the  actual  values  of  the  coordinates,  n  and  5,  are  irrelevant. 
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in  the  same  way  that  r  and  6  may  be  expressed  In  different  units  in  cylin¬ 
drical  coordinates. 

Mow  that  the  values  of  the  coordinates,  n  and  £,  have  been  completely 
specified  on  all  the  boundaries  of  a  closed  field,  it  remains  to  define 
the  values  in  the  interior  of  the  field  in  terms  of  these  boundary  values. 
Such  a  task  immediately  calls  to  mind  elliptic  partial  differential 
equations,  since  the  solution  of  such  an  equation  is  completely  defined 
in  the  interior  of  a  region  by  its  values  on  the  boundary  of  the  region. 
Thus  if  the  coordinates  5  and  n  are  taken  as  the  solutions  of  any  two 
elliptic  partial  differential  equations,  say  L(0  =  0,  D(n)  **  0,  where 
L  and  D  represent  elliptic  operators,  then  €  and  n  will  be  determined 
at  each  point  in  the  interior  of  the  field  by  the  specified  values  on 
the  boundary.  One  condition  must  be  put  on  the  elliptic  system  chosen, 
since  the  same  pair  of  values  (€, h)  must  not  occur  at  more  than  one  point 
in  the  field  or  the  coordinate  system  will  be  ambiguous.  This  condition 
can  be  met  by  choosing  elliptic  partial  differential  equations  exhibiting 
extremum  principles  that  preclude  the  occurrence  of  extrema  in  the  in¬ 
terior  of  the  field. 

This  may  be  illustrated  with  resort  to  the  governing  equation  for 
a  stretched  membrane.  Consider  a  membrane  attached  to  a  flat  plate 
around  a  closed  circuit  of  arbitrary  shape  as  shown  in  Fig.  3.  Now  let 
a  cylinder  of  arbitrary  flat  cross  section  be  pushed  up  through  the  plate, 
stretching  the  membrane  upward.  The  vertical  displacement,  h,  of  the 
membrane  will  be  described  by  Laplace's  equation,  V2h  "  0,  with  h  ■  h^ 
and  h^,  respectively,  on  the  circuits  of  contact  with  the  plate  and  cyl¬ 
inder.  If  equally  spaced  grid  lines  encircling  the  cylinder  had  been 


sioe  view 


TOP  VIEW 


Figure  3.  Illustration  of  extremum  principle 
for  Laplace's  equation 
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drawn  on  the  membrane  before  displacement,  these  lines  would  appear  to 
move  closer  to  the  cylinder  when  viewed  from  above  after  displacement 
of  the  membrane.  None  of  these  lines  would  cross,  however. 

Now  let  pressure  be  applied  on  the  upper  side  of  the  membrane  as 
diagrammed  in  Fig.  4a.  This  will  cause  the  slope  at  the  cylinder  to 
steepen,  with  the  effect  that  the  lines  will  appear  to  be  drawn  even 
closer  to  the  cylinder  but  still  without  crossing.  This  situation  cor¬ 
responds  to  the  Poisson  equation,  V2h  =  p,  where  p  is  the  applied  pressure. 
If  a  variable  pressure  is  applied  on  both  sides  of  the  membrane  to  a 
sufficient  degree,  it  is  possible  to  make  the  membrane  assume  an  S-shape 
as  shown  in  Fig.  4b.  In  this  case  the  encircling  lines  have  crossed, 
and,  consequently,  a  point  on  the  plate  can  no  longer  be  identified  by 
specifying  the  encircling  line  that  it  lies  below  (together  with  a  radial 
ray).  This  latter  case  corresponds  to  a  right-hand  side  of  the  Poisson 
equation  that  is  not  of  one  sign  over  the  entire  membrane,  in  which  case 
the  extremum  principles  of  Poisson's  equation  are  lost. 

Note,  however,  that  if  the  differential  pressure  that  is  applied 
across  the  membrane  is  not  too  large,  the  S-shape  will  not  be  reached. 

In  this  case  the  lines  do  not  cross,  but  rather  the  lines  seem  to  con¬ 
centrate  near  a  line  in  the  interior  of  the  field.  Thus  the  existence 
of  an  extremum  principle  is  a  sufficient  condition  to  prevent  double¬ 
valuedness  in  the  coordinate  system  but  is  not  a  necessary  condition. 

Care  must  be  exercised  in  its  absence,  however. 
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Mathematical  Development 

From  the  discussion  above,  a  logical  choice  of  the  elliptic  gen- 
erating  system  is  Poisson's  equation.  Thus,  based  upon  Fig.  2,  the 
basic  problem  is  to  solve 


£  +  £ 

xx  yy 


-  P 


n  +  n 

xx  yy 


Q 


(1) 


with  boundary  conditions, 


£  *  €^(x,y)  on  ^ 
n  *  constant  *  n  on  T 

(2) 

£  -  £2(x,y)  on  r2 
n  »  constant  =  on  ^ 

The  arbitrary  curve  joining  and  in  the  physical  plane  specifies 
a  branch  cut  for  the  multiple-valued  function,  £(x,y).  Thus  the  values 
of  the  coordinate  functions  x(£,n)  and  y(£,n)  coincide  along  F^ 
and  F^,  and  these  functions  and  their  derivatives  are  continuous  from 

to  Therefore  boundary  conditions  are  neither  required  nor  allowed 

on  and  F  As  previously  noted,  boundaries  with  these  properties 
are  designated  reentrant  boundaries. 

The  functions  P  and  Q  may  be  chosen  to  cause  the  coordinate  lines 
to  concentrate  as  desired,  in  analogy  with  the  membrane  discussed  above. 
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As  discussed  in  Ref.  [7],  negative  values  of  Q  result  in  a  superharmonic 
solution  and  cause  n-lines  to  move  toward  the  n-line  having  the  lowest 
value  of  n,  while  positive  values  have  the  opposite  effect.  Considering 
the  £  solution  to  be  superharmonic  results  in  the  interior  of  the  £  * 
constant  lines  being  rotated  in  a  counterclockwise  direction  in  the  physical 
plane;  whereas  if  the  £-  equation  is  subharmonic,  i.e.,  P  is  positive, 
the  lines  are  rotated  in  the  clockwise  direction.  These  effects 
are  discussed  in  more  detail  below.  It  has  been  found  convenient,  as 
discussed  in  Ref.  [7],  to  redefine  the  control  functions  as 

p  •  7*  <xn2  +  pn2)P 
Q  *  JT  (x?2  +  y^2)Q. 


A  major  purpose  of  this  coordinate  system  control  is  to  concentrate 
lines  in  viscous  boundary  layers  near  solid  surfaces,  and  some  automated 
procedures  for  this  purpose  have  been  developed  (cf.  Ref.  [7]).  Control  is 
also  useful  to  improve  grid  spacing  and  configuration  when  complicated 
geometries  are  involved. 

Since  all  numerical  computations  are  to  be  performed  in  the  rec¬ 
tangular  transformed  plane,  it  is  necessary  to  interchange  the  dependent 
and  independent  variables  in  Eq.  (1).  Using  the  relations  given  in 
Ref.  [3],  Eq.  (1)  becomes 


ax 

ay 


-  2B*Sn  +  ^nn  +  + 

'  2SySn  +  ‘"nr,  +  aPy£  + 


0 

0 


(3) 
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where 


a  m  x2 +  y2 

n  n 

*"Vn+y^n 
Y  *  x|  +  y| 

J  =  Jacobian  of  the  transformation  =  x  v  -  x  y 

V  n  n  C 

with  the  transformed  boundary  conditions 

x  -  r^)  on  r* 

y  =  on  r* 

x  “  f2^’n2^  °n  P2 

y  *  g2(e,n2)  on  r2 

Again  considering  Fig.  2,  the  boundary  functions  f^  f2>  g.^  and  g2 
are  specified  by  the  known  shape  of  the  contours  and  r2  and  the  speci¬ 
fied  distribution  of  K  thereon.  Boundary  data  are  neither  required  nor 
allowed  along  the  reentrant  boundaries  and  Although  the  new 

system  of  equations  is  more  complex  than  the  original  system,  the  boundary 
conditions  are  specified  on  straight  boundaries  and  the  coordinate  spacing 
in  the  transformed  plane  is  uniform.  Computationally,  these  advantages 
far  outweigh  any  disadvantages  resulting  from  the  extra  complexity  of 
the  equations  to  be  solved. 


18 


'nmwi  '  imnntiiiniri  • 


The  boundary-fitted  coordinate  system  so  generated  has  a  constant 


n-  line  coincident  with  each  boundary  in  the  physical  plane.  The  £- 
lines  may  be  spaced  in  any  manner  desired  around  the  boundaries  by 
specification  of  x,y  at  the  equispaced  £-  points  on  the  r£  and  r£ 
lines  of  the  transformed  plane.  As  noted  above,  the  entire  side  boundaries 
are  reentrant  boundaries,  and  thus  neither  require  nor  allow  specification 
of  x,y  thereon. 

Now  the  rectangular  transformed  grid  is  set  up  to  be  the  size 
desired  for  a  particular  problem.  Since  the  values  of  £  and  n  are 
meaningless  in  the  transformed  plane,  the  q-lines  are  assumed  to  run 
from  1  to  the  number  of  n~Hnes  desired  in  the  physical  plane.  Likewise, 
the  £-llnes  are  numbered  1  to  the  number  specified  on  the  boundaries  of 
the  physical  plane.  The  grid  spacing  in  both  the  £  and  q  directions  of 
the  transformed  plane  is  taken  as  unity.  Second-order  central  difference 
expressions  are  used  to  approximate  all  derivatives. 

Only  one  of  a  pair  of  reentrant  boundaries  is  considered  as  a  com¬ 
putation  line  since  the  (x,y)  are  equal  on  both.  As  an  example  of  how 
a  reentrant  boundary  is  handled,  consider  the  grid  in  Fig.  5  where  "o" 
indicates  a  computation  point  and  "A"  a  boundary  point.  The  derivative 
of  x  with  respect  to  £  along  i  =  1  would  be  written  as 


(x2,j  ~  XIMAX-l,j)/2 


(4) 


Again,  it  should  be  stressed  that  all  computations  are  performed 
on  the  rectangular  field  with  square  mesh  in  the  transformed  plane.  The 
resulting  set  of  nonlinear  difference  equations,  two  for  each  point,  is 
solved  by  accelerated  Gauss-Seidel  (SOR)  iteration  using  overrelaxation. 


Figure  5.  Computational  grid  in  transformed  plane 
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Some  discussion  of  this  technique  is  presented  in  Ref.  [8], 

It  might  be  noted  that  both  orthogonal  and  conformal  transformations 
are  special  cases  of  the  generation  of  boundary-fitted  coordinate  systems 
as  the  solutions  of  elliptic  partial  differential  systems.  In  both  of 
these  cases  the  curvilinear  coordinates  satisfy  Laplace's  equation  with 
one  coordinate  constant  on  each  boundary,  and  the  normal  derivative  of 
the  other  coordinate  equal  to  zero  on  each  boundary.  A  conformal  system 
also  requires  a  certain  relation  between  the  range  of  the  two  curvilinear 
coord inates. 

The  same  procedure  may  be  extended  to  regions  that  are  more  than 
doubly  connected,  i.e.,  have  more  than  two  closed  boundaries,  or  equiv¬ 
alently,  more  than  one  body  within  a  single  outer  body.  A  river  reach 
containing  more  than  one  island  would  be  an  example.  One  such  trans¬ 
formation  for  such  a  problem  is  illustrated  in  Fig.  6. 

Types  of  Boundary-Fitted  Coordinate  Systems 

The  above  dLscussion  of  the  generation  of  boundary-fitted  coordinates 
has  centered  around  the  idea  of  using  branch  cuts  to  reduce  multiply 
connected  regions  to  simply  connected  ones  in  the  transformed  plane. 

An  example  using  branch  cuts  is  sketched  in  Fig.  7 .  Here  the  body  in 
the  field  transforms  to  the  entire  bottom  boundary  of  the  transformed 
plane,  while  the  entire  surrounding  boundary,  1-2-3-4-5-6, 
transforms  to  the  top  boundary  of  the  transformed  plane.  The  sides  of 
the  transformed  plane  are  reentrant  boundaries,  corresponding  to  the  cut, 
8-1  and  7-6,  in  the  physical  field.  Thus,  in  the  difference  equations 
points  lying  just  to  the  right  of  the  right  boundary  are  identical  with 
corresponding  points  just  to  the  right  of  the  left  boundary.  This  is 
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TRANSFORMED  PLANE 

Figure  6.  Boundary-fitted  coordinates  for  a  river 
containing  two  islands 


22 


TRANSFORMED  PLANE 


Figure  7.  Example  of  coordinates  generated  using  a 
branch  cut.  Placement  of  body  is  such  that  sides 
are  reentrant  boundaries 
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the  same  type  of  circumstance  that  occurs  with  the  familiar  cylindrical 
coordinate  system,  where  6  =  361°  is  the  same  point  as  6  9  1°.  Similarly, 
points  just  outside  the  left  boundary  are  coincident  with  points  just 
inside  the  right  boundary. 

Many  variations  of  this  type  of  coordinate  system  can  be  produced, 
cf.  Ref.  [3],  For  instance,  the  transformed  plane  corresponding  to  the 
same  physical  field  shown  in  Fig.  7  can  be  rearranged  as  shown  in  Fig. 

8.  Now  the  reentrant  boundary,  corresponding  to  the  cut,  is  located  on 
a  portion  of  the  bottom  of  the  transformed  plane.  The  coordinate  lines 
that  result  from  these  two  types  of  arrangements  of  the  transformed  plane 
are  shown  on  each  of  the  figures.  As  with  all  the  boundary-fitted  coor¬ 
dinate  systems,  the  grid  is  square  in  the  transformed  plane  regardless 
of  the  line  configuration  in  the  physical  plane. 

Multiple-body  fields  can  also  be  transformed  to  simply  connected 
regions,  an  example  of  which  is  shown  in  Fig.  9  .  Again  there  are  many 
different  possible  arrangements  of  the  transformed  plane,  all  of  which 
are  created  by  sliding  the  boundary  segments  around  the  rectangular 
boundary  of  the  transformed  plane.  A  number  of  examples  are  given  in 
kef.  [3]  and  Ref.  [8]. 

The  other  type  of  coordinate  system  transformation  available  leaves 
the  multiplicity  of  the  region  unchanged.  In  this  case,  bodies  in  the 
interior  of  the  physical  field  are  transformed  to  rectangular  slabs  or 
even  slits  in  the  transformed  plane.  Three  different  possibilities  are 
shown  in  Fig.  10  for  the  physical  plane  shown  in  Fig.  7,  In  the  case  of 
slits,  the  physical  coordinates  and  solution  variables  in  general  have 
different  values  at  points  on  the  two  sides  of  the  slit,  even  though  such 
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PHYSICAL  PLANE 


Figure  8.  Example  of  coordinates  generated  using  a 
branch  cut.  Placement  of  body  is  such  that  reentrant 
boundaries  lie  on  bottom  line  of  the  transformed 

plane 
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-REENTRANT  BOUNDARIES 


TRANSFORMED  PLANE 


Figure  9.  Coordinates  generated  for 
multiple-body  field 


points  are  coincident  in  the  transformed  plane.  This  does  not  introduce 
any  approximations,  but  simply  adds  a  little  more  bookkeeping  to  the 
code.  Fields  with  more  than  one  body  in  the  interior  simply  result  in 
a  like  number  of  slabs  and/or  slits  in  the  transformed  plane. 

Comparison  of  all  of  the  above  figures  shows  that  different  types 
of  transformation  may  be  more  appropriate  for  different  physical  config¬ 
urations.  A  further  example  of  this  is  the  configuration  in  Fig.  11, 
shown  with  three  variations.  Generally,  the  slit/slab  form  is  more 
appropriate  for  channel-like  physical  configurations  having  bodies  in 
the  interior,  while  the  other  form  works  particularly  well  for  "unbounded" 
regions  involving  external  flow  about  bodies  and  for  regions  having  an 
outer  boundary  that  forms  a  continuous  circuit  without  pronounced  corners 
around  the  field.  The  slab  is  generally  superior  to  the  slit  unless 
the  boundary  has  a  sharp  point.  The  case  of  a  single  channel  without 
any  interior  bodies  is  the  same  in  either  form.  An  example  of  a  river 
reach  containing  two  islands,  using  horizontal  slits  rather  than  the 
branch  cuts  previously  presented  in  Fig.  6, is  given  in  Fig.  12. 

Data  Required  for  Generation  of  Boundary-Fitted  Coordinates 

The  basic  input  or  data  required  to  generate  a  boundary-fitted 
coordinate  system  are  the  physical  coordinates  of  points  on  the  boundaries. 
For  example,  with  reference  to  Fig.  7,  the  coordinates  of  points  on  the 
body  from  8  around  to  7  would  be  required,  with  these  points  being 
spaced  in  any  manner  desired  as  long  as  there  is  a  continuous  progression 
from  8  to  7.  Similarly,  the  (x,y)  values  for  points  on  the  outer  boundary 
from  1  to  2,  etc.,  on  around  to  6  would  be  required.  Again  these  points 
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TRANSFORMED  PLANE 


Figure  11.  Comparison  of  coordinates  without  interior 
body  and  with  slit/slab  interior  body 


PHYSICAL  PLANE 


TRANSFORMED  PLANE 


Figure  12.  Coordinates  generated  with  slits  for 
river  with  two  islands 


may  be  spaced  around  the  boundary  as  desired,  with  no  restriction  as  to 
how  many  points  lie  on  each  boundary  segment,  e.g.,  between  1  and  2  or 
between  4  and  5,  provided  that  only  the  total  number  of  points  from  1 
around  to  6  is  the  same  as  from  8  to  7.  The  coordinates  of  points  must 
be  specified  on  the  entirety  of  these  lines.  The  coordinates  of  points 
on  reentrant  segments  of  the  boundary  in  the  transformed  plane,  e.g.,  1 
to  8  and  6  to  7,  are  not  specified  but  are  free  to  be  determined  by  the 
solution. 

Similarly,  with  reference  to  Fig.  10a,  the  coordinates  of  outer 
boundary  points  are  required  In  the  slab/slit  transformations.  In 
addition,  body  points  from  6  to  1  on  the  lower  half  of  the  body  and 
from  1  to  6  on  the  top  half  are  required.  No  calculations  would  be 
made  on  the  slab  sides  of  Figure  10c  or  slits  of  Figures  10a  and  10b 
since  values  at  such  points  are  fixed.  Points  in  the  interior  of  a 
slab  are  irrelevant.  As  always,  points  may  be  spaced  as  desired  around 
the  bodies  and  outer  boundary  segments. 

Computer  Time  Required  for  Generation  of  Boundary-Fitted  Coordinates 

Ref.  [  8 ]  Indicates  that  the  typical  time  required  to  generate  a 
one-body  coordinate  system  without  coordinate  system  control  (the 
functions  P  and  Q  are  set  to  zero)  is  about  2  min  on  a  UNIVAC  1106  com¬ 
puter  for  a  70  x  30  field  (70  points  on  the  body).  If  P  and  Q  are  not 
zero,  so  that  the  spacing  of  coordinate  lines  is  controlled,  the  computation 
time  increases.  Multiple-body  coordinate  systems  typically  require  about 
6  min  for  a  70  x  40  field.  If  these  same  computations  were  to  be  made 
on  a  CDC-7600  computer,  the  times  quoted  above  would  be  reduced  by  perhaps 
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an  order  of  magnitude  or  more.  Therefore,  the  cost  of  generating 
boundary-fitted  coordinate  systems  for  use  in  numerical  models  will 
be  generally  insignificant. 

COORDINATE  SYSTEM  CONTROL 

Control  of  the  coordinate  line  spacing  in  the  field  can  be  exercised 
through  the  non-zero  values  given  to  the  Laplacian  of  the  curvilinear 
coordinates  as  in  Eq.  (1),  as  noted  above.  With  a  zero  Laplacian,  the 
lines  tend  to  be  closely  spaced  near  convex  segments  and  more  widely 
spaced  near  concave  segments.  A  negative  value  of  the  Laplacian  causes 
the  lines  to  move  toward  lower  values  of  the  curvilinear  coordinate. 

Attraction  to  Other  Coordinate  Lines  and/or  Points 

This  effect  is  utilized  as  in  Ref.  [  8 ]  to  achieve  attraction  of 
coordinate  lines  to  other  coordinate  lines  and/or  points  by  taking  the 
form  of  the  control  functions  to  be 

n 

p(£.n)  *  -  l  a±  sign(C  -  ^i)exp(-ci  |  £  -  cj) 
i-1 

(5) 

-  I  sign(C  -  Ci)exp{-d1[(C  -  q)2  +  (n  -  ^ ) 2 1  } 
i-1 

and  an  analogous  form  for  ()(£,q)  with  £  and  q  interchanged.  The  effects 
of  such  control  is  illustrated  in  Refs.  [  7]  and  [8],  The  efficacy 
of  control  to  improve  the  accuracy  of  a  physical  solution  done  on  the 
coordinate  system  has  been  noted. 

In  the  P  function,  the  effect  of  the  amplitude,  a^,  is  to  attract 
C  -  coordinate  lines  toward  the  5^-line,  while  the  effect  of  the  amplitude 
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is  to  attract  5-lines  toward  the  single  point  Note  that 

this  attraction  to  a  point  is  actually  attraction  ot  £-lines  to  a  point 
on  another  5-line,  and,  as  such,  acts  normal  to  the  £-line  through  the 
point.  There  is  no  attraction  of  q-lines  to  this  point  via  the  P 
function.  In  each  case  the  range  of  the  attraction  effect  is  determined 
by  the  decay  factors,  c^  and  d^.  With  the  inclusion  of  the  sign  changing 
function,  the  attraction  occurs  on  both  sides  of  the  £-line,  or  the 
(£i»  Of)  point,  as  the  case  may  be.  Without  this  function,  attraction 
occurs  only  on  the  side  toward  increasing  5.,  with  repulsion  occurring  on 
the  other  side.  A  negative  amplitude  simply  reverses  all  of  the  above- 
described  effects,  i.e.,  attraction  becomes  repulsion  and  vice  versa. 

The  effect  of  the  Q.  function  of  n~lines  follows  analogously.  It  should 
be  noted  that  P  and  Q.  are  discontinuous  because  of  the  sign  function  and 
are  equal  to  sums  of  second  derivatives.  As  a  consequence,  the  coordinates 
have  continuous  first  derivatives  but  discontinuous  second  derivatives 
at  controlled  locations. 

In  the  case  of  a  boundary  that  is  an  n-line,  positive  amplitudes 
in  the  Q  function  will  cause  n-lines  off  the  boundary  to  move  closer 
to  the  boundary,  assuming  that  n  increases  off  the  boundary.  The  effect 
of  the  P  function  will  be  to  alter  the  angle  at  which  the  5-lines  inter¬ 
sect  the  boundary,  since  the  points  on  the  boundary  are  fixed,  with  the 
5-lines  tending  to  lean  in  the  direction  of  decreasing  5.  If  the  boundary 
is  such  that  n  decreases  off  the  boundary,  then  the  amplitudes  in  the  Q 
function  must  be  negative  to  achieve  attraction  to  the  boundary.  In 
any  case,  the  amplitudes  a^  cause  the  effects  to  occur  all  along  the 
boundary,  while  the  effects  of  the  amplitudes  b^  occur  only  near  se¬ 
lected  points  on  the  boundary. 
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Attraction  to  Space  Curves  and/or  Points 


If  the  attraction  line  and/or  the  attraction  points  are  in  the 
field,  rather  than  on  a  boundary,  then  the  attraction  is  not  to  a  fixed 
line  or  point  in  space,  since  the  attraction  line  or  points  are  themselves 
solutions  of  the  system  of  equations,  the  functions  P  and  Q_  being  functions 
of  the  variables  £  and  q.  It  is,  of  course,  also  possible  to  take  these 
control  functions  as  functions  of  x  and  y,  instead  of  £  and  q ,  and  achieve 
attraction  to  fixed  lines  and/or  points  in  the  physical  field.  This 
case  becomes  somewhat  more  complicated,  since  it  must  be  ensured  that 
coordinate  lines  are  not  attracted  parallel  to  themselves.  The  following 
development  was  given  in  Ref.  [9]. 

Recall  that  in  the  above  discussion,  q-lines  are  attracted  to  other 
q- lines  ,  and  £- lines  are  attracted  to  other  £- lines  .  It  is  unreasonable, 
of  course,  to  attempt  to  attract  q-  lines  to  £- lines  ,  since  that  would 
have  the  effect  of  collapsing  the  coordinate  system: 


When,  however,  the  attraction  is  to  be  to  certain  fixed  lines  in 
x-y  space,  defined  by  curves  y  ■  f(x),  care  must  be  exercised  to  avoid 
attempting  to  attract  q-  or  £-lines  to  specified  curves  that  cut  the 
n-  or  £-lines  at  large  angles.  Thus,  in  the  figure  below, 
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5 -line 


it  is  unreasonable  to  attract  £-lines  to  the  curve  f(x),  while  it  is 
natural  to  attract  the  n-lines  to  f(x). 

However,  in  the  general  situation,  the  specified  line  f(x)  will  not 
necessarily  be  aligned  with  either  a  £-  or  n-line  along  its  entire  length. 
Since  it  is  unreasonable  to  attract  a  line  tangentially  to  itself,  some 
provision  is  necessary  to  decrease  the  attraction  to  zero  as  the  angle 
between  the  coordinate  line  and  the  given  line  f(x)  goes  to  90°.  This 
can  be  accomplished  by  multiplying  the  attraction  function  by  the  cosine 
of  the  angle  between  the  coordinate  line  and  the  line  f(x).  It  is  also 
necessary  to  change  the  sign  on  the  attraction  function  on  either  side 
of  the  line  f(x).  This  can  be  done  by  multiplying  by  the  sine  of  the 
angle  between  the  line  f(x)  and  the  vector  to  the  point  on  coordinate 
line. 

These  two  purposes  can  be  accomplished  as  follows.  Let  a  general 
point  on  the  £-line  be  located  by  the  vector  R(x,y),  and  let  the  attrac¬ 
tion  line  y  •  f(x)  be  specified  by  the  collection  of  points  SCx^.y^), 
i  *  1,  2,  n.  Let  the  unit  tangent  to  the  attraction  line  be 
t(*  fy  ),  and  the  unit  tangent  to  a  £-line  be  t^. 


The  control  functions  P(x,y)  and  0(x, y)  may  then  be  logically  taken  as 


n 

p(x,y)  =  -  I  ai(ti 

i=l 


(e\  U.x  (R  -  S3]  •  k 
I  }  - |R  “Sj - ^  expC-djR  -  S.|) 


(6) 


°S*,y)  =  -  [  ajL(t1 
i=l 


[t±x  (R  -Sj)] 

|R -'sj 


k 

—  exp(-di|R  -  Si|) 


where  k  is  the  unit  vector  normal  to  the  two-dimensional  plane.  These 
relations  are  evident  from  the  figure  below: 


angle  between  the  C-line  and  the  attraction  line  approaches  90°.  The 
cross  product  term  changes  the  sign  of  the  control  function  on  either 
side  of  the  attraction  line  to  produce  attraction  on  both  sides  of  the 


line.  Again  the  strength  and  range  of  the  attraction  are  determined  by 
the  amplitude,  a^,  and  the  decay  factor,  d^  respectively. 
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These  functions  depend  on  x  and  y  through  both  R  and  or  and 

thus  must  be  recalculated  at  each  point  as  the  iterative  solution  proceeds. 
This  form  of  coordinate  control  will  therefore  be  more  expensive  than 
that  based  on  attraction  to  other  coordinate  lines. 

There  is  no  real  distinction  between  "line"  and  "point"  attraction 
with  this  type  of  attraction.  "Line"  attraction  here  is  simply  attraction 
to  a  group  of  points  that  form  a  line  f(x).  If  line  attraction  is  speci¬ 
fied,  then  the  tangent  to  the  line  f(x)  is  computed  from  the  adjacent 
points  on  the  line.  If  point  attraction  is  specified,  then  the  "tangent" 
must  be  input  for  each  point.  The  tangents  to  the  coordinate  lines  are 
computed  from  the  relations  given  in  Ref.  [3]. 

Control  Functions  from  Boundary-Po int  Distributions 

With  the  Laplacians  of  the  coordinates  equal  to  zero,  the  line 
spacing  in  the  field  will  not  be  greatly  affected  by  the  distribution 
of  the  boundary  points,  except  very  near  the  boundaries.  In  fact,  if 
the  control  functions  are  not  consistent  with  the  boundary  point  dis¬ 
tribution,  very  large  changes  in  the  metric  coefficients  will  occur  near 
the  boundaries.  Values  of  the  control  functions  may  be  determined  from 
the  ID  boundary  point  distribution  such  that  the  line  spacing  in  the 
field  will  generally  follow  that  on  the  boundary.  This  concept  was  in¬ 
troduced  in  Ref.  [10  ]  and  is  discussed  in  Ref.  [7]  as  generalized  to  3D 
in  Ref.  [  11  ] .  However,  in  the  use  of  control  functions  that  are 
ID,  it  should  be  noted  that  excessive  concentration  of  lines  can  occur 
near  sharp  convex  corners  as  discussed  in  Ref.  [7]. 

With  Eq.  (3)  evaluated  in  ID  on  a  straight  n-line  conincident  with 

the  x-axis,  we  have,  since  x  =  y  *  0  in  this  case, 

n 
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ax^  --aP(£)x£  (7) 

The  reason  for  the  choice  of  the  form  of  the  control  functions  in  Eq.  3 
becomes  clear,  since  a  cancels  from  this  equation  to  leave 

P(£)  =  -x^/x^  (8) 

Thus  the  control  function,  P(£)»  can  be  determined  from  the  specified 
boundary  point  distribution,  x(£).  Generalizing,  x  is  replaced  by  arc 
length  along  the  £-line  ,  and  the  effect  will  be  qualitatively  the  same 
when  this  line  is  curved,  (Compare  Ref.  [7]  for  more  detail.) 

If  this  value  of  the  control  function  is  then  used  throughout  the 
field,  the  £-line  distribtuion  in  the  field  will  generally  follow  the 
specified  distribution  of  the  end  points  of  these  lines  on  the  boundary. 
With  different  point  distributions  on  two  boundaries,  values  of  the 
control  function  P(£»n)  in  the  field  between  can  be  determined  by  ID 
interpolation  in  n  between  the  values  determined  in  the  above  manner  on 
the  two  q-line  boundaries.  An  analogous  development  applies  for  the 
determination  of  the  control  function  (£(£,n)  from  interpolation  in 
between  ID  evaluations  on  two  £-line  boundaries.  This  interpolation 
was  introduced  in  Ref.  [  12  ]  in  a  2D  coordinate  system. 

SYSTEM  CONFIGURATION 

In  the  present  model,  the  physical  field  may  have  both  external  and 
internal  boundaries  of  arbitrary  shape.  The  field  in  the  transformed 
plane  is  rectangular  with  rectangular  holes  corresponding  to  any  internal 
boundaries.  This  configuration  is  illustrated  in  Fig.  13.  Boundary 
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PHYSICAL  PLANE 


Figure  13.  Example  of  coordinates  generated  in  a 
field  containing  a  jetty  and  an  island 


I 


Intrusions  may  be  transformed  either  to  portions  of  the  rectangular  outer 
boundary  of  the  transformed  region,  as  in  Fig.  13,  or  to  slabs  protruding 
inward  from  this  boundary  as  in  Fig.  14.  A  general  discussion  of  possible 
configurations  is  given  in  Ref.  [3].  Various  outlet  shapes  and  locations, 
as  well  as  internal  obstacles  and  boundary  protrusions  such  as  weirs, 
can  be  treated  by  the  same  code  with  only  changes  in  the  input.  This 
input  consists  of  the  physical  cartesian  coordinates  of  the  points  se¬ 
lected  on  each  segment  of  the  physical  boundaries.  A  small  front-end 
code  was  written  to  provide  certain  line  segments  (linear,  quadratic, 
and  cubic  polynomials)  with  linear  or  exponential  distributions  thereon 
automatically. 

The  code  automatically  calculates  control  functions  P(£,n)  and 
<2(£,n)  for  the  coordinate  generation  equations  (3)  from  the  boundary 
point  distribution  as  discussed  above.  These  functions  are  calculated 
from  the  ID  relations  on  each  boundary  segment  and  are  interpolated 
linearly  into  the  field  between  opposing  boundary  sections  in  the 
transformed  plane. 

In  addition,  attraction  of  coordinate  lines  to  other  coordinate 
lines  and/or  points,  and  to  specified  lines  and/or  points  in  space,  also 
discussed  above,  is  provided  through  input  quantities.  This  input 
consists  of  the  coordinate  lines  and/or  points,  and  the  specified  space 
curves  and/or  points,  to  which  the  attraction  is  to  be  made  and  the  ampli¬ 
tudes  and  decay  factors  for  the  corresponding  attractions. 
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Several  examples  of  coordinate  systems  produced  by  this  code  are 
given  in  Figs.  15-19.  Examples  of  applications  of  such  systems  appear 
in  Ref.  [4]-[6],  Two  further  examples,  together  with  complete  input 
listings  for  the  code,  follow  the  description  of  the  code  in  Part  B. 
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PART  B 


COORDINATE  CODE 

The  present  code  (WESC0R)  differs  from  a  previous  version  (TOMCAT) 
described  in  Ref.  [8]  in  that  the  latter  does  not  provide  for  slits  and/or 
slabs  in  the  interior  of  the  transformed  plane.  Also,  branch  cuts  (if  used) 
in  the  present  code  are  restricted  to  the  entire  left  and  right  sides  of 
the  outer  rectangle  in  the  transformed  region.  Finally,  the  present  code 
includes  a  more  extensive  means  of  coordinate  line  control,  involving 
attraction  to  space  lines/or  points  and  also  involving  determination  from 
boundary  point  distributions. 

The  code  for  the  numerical  generation  of  the  boundary-fitted  coor¬ 
dinate  system  from  the  equations  of  Part  A,  together  with  a  front-end 
code  to  generate  boundary  point  distributions  and  a  plot  code,  is  discussed 
below.  These  codes  were  implemented  on  the  CRAY-1  computer  at  the  Air 
Force  Weapons  Laboratory,  Kirtland  AFB,  New  Mexico. 

WESC0R  (Coordinate  System) 

This  code  generates  the  boundary- f itted  coordinate  system  by  solving 
a  set  of  elliptic  partial  differential  equations  by  SOR  iteration  as 
discussed  in  Part  A.  Attraction  of  coordinate  lines  to  other  coor¬ 
dinate  lines  and/or  points, and  to  specified  lines  and/or  points  in  space, 
is  included.  The  shape  and  configuration  of  the  boundary  are  arbitrary, 
except  that  the  outer  boundary  must  be  closed.  There  may  be  an  arbitrary 
number  of  internal  closed  boundaries  transforming  to  either  slits  or 
slabs  as  discussed  in  Part  A. 

The  input  to  this  code  consists  of  the  point  distribution  on  the 
boundary  of  the  region,  several  quantities  in  connection  with  the  control 
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of  the  coordinate  line  spacing,  and  the  parameters  associated  with  the 
iterative  solution  process.  This  input  is  described  in  detail  below. 

The  file  output  from  the  code  LINES  can  be  used  directly  as  a  part  of 
the  input  to  this  code  from  file  10.  A  simplified  flow  chart  of  WESC0R 
is  shown  in  Figure  20. 

Boundary  Configurations 

Arrays .  The  dependent  variable  field  arrays  are  X  and  Y,  which  contain 
the  cartesian  coordinates  (x,y)  for  each  grid  point.  The  indices  (I,J) 
of  these  arrays  correspond  to  the  curvilinear  coordinates  (  £,  q) .  and 
run  from  1  to  IMAX  and  JMAX,  respectively.  The  increments  A£  and  An 
in  the  difference  expressions  are  thus  equal  to  unity  by  construction. 

(  These  increments  cancel  from  all  the  difference  equations  and  are  thus 
irrelevant . ) 

In  order  to  treat  slit  configurations,  for  which  a  closed  interior 
boundary  in  the  physical  region  is  collapsed  to  a  slit  in  the  transformed 
region,  there  are  four  other  coordinate  arrays,  XL,  YL  and  XU,  YU,  which 
contain  the  cartesian  coordinates  on  the  two  sides  of  the  slit.  The 
first  index  of  these  arrays  corresponds  to  the  location  of  the  point 
relative  to  the  left  end  of  horizontal  slits,  or  relative  to  the  lower 
end  of  vertical  slits,  this  end  index  being  designated  unity.  The  other 
index  Indent  if ies  the  particular  slit.  For  horizontal  slits  the  coor¬ 
dinates  on  the  lower  side  are  in  XL  and  YL,  while  those  on  the  upper 
side  are  in  XU  and  YU.  Vertical  slits  have  the  coordinates  on  the  left 
side  in  XL  and  YL,  and  those  on  the  right  side  in  XU  and  YU. 

There  is  also  a  field  array  LSLIT(I,J)  containing  the  point  type 
for  each  point.  This  array  identifies  each  point  as  being  on  a  slit. 
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adjacent  to  a  silt,  on  a  slab  side,  on  an  outer  boundary,  in  the  field, 
or  out  of  the  computation  region  (inside  a  slab) ,  as  illustrated  on  the 
diagram  below: 
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The  coordinate  system  control  functions  P  and  0  for  each  point 
are  contained  in  the  field  arrays  RXI(I,J)  and  RETA(I,J),  respectively. 
There  are  also  arrays  RXIL,  RETAL  and  RX1U,  RETAU,  analogous  to  the 
array  XL,  etc.,  discussed  above,  which  contain  the  values  of  these 
functions  on  the  two  sides  of  the  slits.  The  acceleration  parameters 
for  the  iteration  at  each  point  are  in  the  field  array  WACC(I,J). 

Configuration  types.  The  cartesian  coordinates  of  the  points  on  the 
entire  boundary  of  the  physical  region,  i.e.,  the  closed  outer  boundary 
and  any  internal  boundaries,  must  be  input.  There  are  two  basic  types 
of  overall  configuration  included  in  the  code.  In  one  the  connectivity 
of  the  transformed  region  is  the  same  as  that  of  the  physical  region, 
i.e.,  the  closed  outer  boundary  of  the  physical  region  corresponds  to  a 
closed  outer  boundary  of  the  transformed  region.  With  the  other  type, 
one  branch  cut  is  introduced  in  the  physical  region  so  that  the  closed 
outer  boundary  and  one  inner  boundary  of  the  physical  region  transform 
to  the  bottom  and  top  of  a  rectangle  forming  the  outer  boundary  of  the 
transformed  region.  The  left  and  right  sides  of  the  transformed  region 
then  correspond  to  the  branch  cut  in  the  physical  region.  Points  on 
these  sides  therefore  are  not  input  but  rather  are  calculated  as  part  of 
the  solution. 

Rectangular  outer  boundary.  If  the  outer  boundary  of  the  physical 
region  is  to  correspond  to  a  rectangle  forming  the  outer  boundary  of  the 
transformed  region,  then  the  points  on  this  boundary  can  be  input  in 
clockwise  succession  around  the  outer  rectangle  of  the  transformed  region 
as  in  the  diagram  below.  If  the  outer  boundary  of  the  physical  region 
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is  a  circle,  then  the  points  on  this  circle  can  be  generated  internally 
by  the  code,  requiring  input  only  of  the  radius  (YINFIN)  and  cartesian 
coordinates  of  the  center  (X0INF,Y0INF)  of  the  circle,  together  with 
the  cartesian  coordinates  of  the  angular  position  (AINFIN)  and  indices 
(1NFXI,1NFETA)  of  the  point  at  which  the  clockwise  succession  of  points 
around  the  outer  rectangle  is  to  start,  and ,/the  total  number  of  points 
on  the  circle  (NINF) .  As  above,  the  points  will  be  placed  in  clockwise 
succession  around  the  circle  or  boundary  of  the  physical  region  and 
the  rectangular  boundary  of  the  transformed  region.  The  treatment  of 
the  outer  boundary  is  determined  by  the  input  parameter  IBNDRY. 

An  alternative  procedure  for  inputting  the  outer  boundary  is  to  input 
each  straight  segment  of  this  boundary  of  the  transformed  region  as  a 
slab  side  in  the  manner  described  below  for  internal  boundaries. 

Internal  boundaries  (slits/slabs).  Internal  boundaries  in  overall 
configurations  of  the  former  type  introduced  above  correspond  to  either 
slits  or  slabs  in  the  transformed  region: 
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In  the  case  of  slits,  the  points  are  Input  In  clockwise  succession 
beginning  at  the  right  end  for  horizontal  slits  or  counter-clockwise  be¬ 
ginning  with  the  top  for  vertical  slits,  and  are  placed  in  the  arrays  XL, 
etc.,  described  above.  For  slabs,  the  four  sides  are  input  independently 
and  the  succession  of  points  may  be  in  either  direction  on  each  side. 

In  fact,  it  is  not  even  necessary  for  the  four  sides  of  one  slab  to  be  input 
in  succession;  the  sides  of  all  slabs  in  the  field  may  be  placed  in  any 
order  in  the  input.  The  coordinates  of  the  points  on  slab  sides  are 
placed  directly  in  the  field  arrays  X  and  Y.  This  input  of  boundary 
segments  corresponding  to  slits  or  slabs  is  accomplished  as  follows. 

For  horizontal  slits,  the  5-indices(I)  0f  the  left  and  right  ends 
are  placed  in  the  arrays  LBl  and  LB2,  respectively.  The  n-index  (J)  of 
the  entire  slit  or  slab  side  is  placed  in  the  array  LB3.  In  the  case 
of  vertical  slits,  the  n- indices  (J)  of  the  bottom  and  top  go  in  LBl 
and  LB2,  while  the  5-index  (I)  goes  in  LB3.  Slab  sides  are  treated  in 
the  same  manner  except  that,  since  the  points  thereon  may  be  input  in 
either  direction,  LBl  and  LB2  contain  the  indices  of  the  end  points  of 
the  side  in  either  order,  i.e.,  LBl  may  exceed  LB2.  The  points  are  input 
from  LBl  to  LB2. 

For  both  slits  and  slab  sides,  a  flag  is  placed  in  an  array  LTYPE 
to  designate  the  segment  as  a  slit  or  slab  side  in  horizontal  or  vertical 
orientation: 

+1  horizontal  slit 
+2  vert ical  si it 
-1  horizontal  slab  side 
-2  vertical  slab  side 

The  code  computes  the  number  of  points  on  the  slit  or  slab  side  from  the 
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values  of  LB1  and  LB2  and  places  this  value  In  the  array  LPT.  All  of  these 
arrays  are  single-dimension  arrays,  there  being  one  set  of  parameters 
for  each  slit  or  slab  side.  The  total  number  of  slits  and  slab  sides. 


including  those  on  the  outer  boundary  as  described  below,  is  specified 
by  the  input  parameter  NBDY. 

Outer  boundary  Intrusions.  As  noted  above,  the  outer  boundary  can  be 
input  In  segments  as  slab  sides.  This  is  illustrated  below. 


This  is  done  just  as  described  above  for  internal  boundaries  except  that 
values  of  -11  and  -12,  respectively,  are  input  for  LTYPE  for  horizontal 
and  vertical  segments  of  the  outer  boundary. 

Branch  cut.  With  the  other  type  of  overall  configuration,  involving 
a  branch  cut,  the  outer  boundary  and  the  internal  boundary  connected  to 
the  cut  are  both  input  clockwise  from  the  points  joined  by  the  cut.  As 
noted  above,  these  points  are  placed  on  the  top  and  bottom  of  the  rec¬ 
tangle  forming  the  outer  boundary  of  the  transformed  region.  This  type 
of  configuration  is  elected  through  the  input  parameter  NREN.  Additional 
Internal  boundaries  can  be  input  as  either  slits  or  slabs  exactly  as 
described  above. 

Boundary  input.  Provision  is  made  for  reading  the  boundary  points 
either  from  card  images  (x  and  y  for  one  point  to  a  card  in  2F10.0 
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format)  or  from  the  output  of  the  LINES  code  described  below,  as  de¬ 
termined  by  the  input  parameter  ISLIT.  The  outer  boundary  must  be  input 
as  segments  of  slab  sides  if  this  boundary  is  included  on  the  output  of 
the  LINES  code. 

Control  Functions 

Coordinate  system  control  is  included  through  both  the  attraction 
of  coordinate  lines  to  other  coordinate  lines  and/or  points  and  to  speci¬ 
fied  lines  and/or  points  in  the  physical  region,  as  described  in  Part  A. 
(For  completeness,  provision  is  made  for  repulsion  as  well  as  attraction.) 

Attraction  to  coordinate  lines  and/or  points.  The  first  of  these 
requires  the  input  of  the  index  (indices)  of  the  curvilinear  coordinate 
line,  together  with  the  associated  attraction  amplitude  and  decay  factor, 
for  each  line  (point)  to  which  the  attraction  is  made.  For  attraction 
to  lines,  the  index,  amplitude,  and  decay  factor  are  placed  in  the  arrays 
JLN,  ALN,  and  DLN,  respectively,  while  for  attraction  to  points,  the 
corresponding  arrays  IPT,  JPT,  APT,  and  DPT  are  used. 

Attraction  to  space  lines  and/or  points.  For  attraction  to  specified 
lines  and/or  points  in  space,  the  input  is  similar  in  regard  to  the  ampli¬ 
tude  and  decay  factors,  using  the  arrays  APT  and  DPT.  It  is  necessary, 
of  course,  to  also  input  the  cartesian  coordinates  of  the  points  on  the 
line,  or  the  Isolated  points,  to  which  the  attraction  is  made.  These 
coordinates  are  placed  in  the  arrays  XPT  and  YPT.  For  attraction  to 
points,  it  is  also  necessary  to  input  the  components  of  a  vector  normal 
to  the  desired  direction  of  the  attraction  for  each  point,  these  com¬ 
ponents  being  placed  in  the  arrays  VEC1  and  VEC2. 
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Effect  of  boundary  point  distribution.  In  addition  to  the  above  types 
of  attraction,  the  control  functions  also  include  the  effect  of  the 
boundary  point  distribution  discussed  in  Part  A.  This  is  done  by  evalu¬ 
ating  one  of  the  control  functions  on  each  boundary  segment  in  the 
transformed  region  (P  on  n-lines,  Q.  on  £-lines)  from  the  one-dimen¬ 
sional  relations  in  terms  of  arc  length  discussed  in  Part  A.  These  values 
are  placed  in  the  arrays  RXI  and  RETA,  except  for  slits  where  the  arrays 
RXIL,  etc.,  are  used  in  the  manner  described  above  for  XL,  etc.  Values 
of  the  control  functions  in  the  field  are  then  interpolated  linearly 
between  facing  boundary  segments,  P  being  interpolated  vertically  and 
Q.  horizontally,  as  illustrated  in  the  following  diagram. 


J - F 


This  evaluation  is  done  first  and  then  the  contributions  to  the  control 
functions  from  the  line  and  point  attraction  is  added  to  the  arrays 
RXI  and  RETA  in  the  field. 

Iterative  Solution 

Initial  guess.  The  initial  guess  for  the  values  of  the  cartesian  co¬ 
ordinates  in  the  field,  i.e.,  the  values  in  the  arrays  X  and  Y  in  the 
field,  that  is  used  to  start  the  iterative  solution  is  obtained  by  the 
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same  type  of  interpolation  between  facing  segments  described  above  for 
the  control  functions,  except  that  both  X  and  Y  are  interpolated  between 
the  pair  of  facing  segments  with  the  smallest  separation  in  the  transformed 
region.  Thus  values  at  point  1  in  the  figure  below  would  be  obtained 
by  horizontal  interpolation,  but  at  2  the  interpolation  would  be  vertical. 


1 

1 

1 

• 

1 

1 

• 

■ 

L_J _ S _ □ 

Since  very  strong  control  functions  can  sometimes  make  the  conver¬ 
gence  of  the  iterative  solution  difficult  in  complicated  configurations, 
provision  is  made  for  first  converging  the  field  with  the  control  functions 
set  to  zero  and  then  re-converging  in  steps  as  these  functions  are 
increased  to  full  value.  Actually  this  feature  is  rarely  needed. 

Acceleration  parameters.  As  discussed  in  Part  A  the  solution  for  the 
cartesian  coordinates  in  the  field  is  done  by  SOR  iteration.  Either  a 
uniform  value  of  the  acceleration  parameter  can  be  input  as  R(l)  or  the 
code  will  calculate  a  locally  optimum  value  at  each  point  in  the  field, 
these  values  being  placed  in  the  field  array  WACC.  This  calculation 
is  discussed  in  Ref.  [8],  where  it  is  noted  that  the  values  obtained  are 
not  truly  optimum  in  all  cases.  Therefore  this  provision  has  not  been 
found  to  be  as  generally  efficient  as  simply  using  a  uniform  valuer  since  the 
calculation  of  the  acceleration  parameter  involves  a  square  root  and 
hence  is  time-consuming.  The  uniform  value  should  be  around  1.85  for 
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large  fields.  This  value  should  be  decreased  for  strong  control  functions 
or  small  fields. 

Iterative  process.  The  iteration  continues  until  either  the  magni¬ 
tude  of  the  changes  in  the  values  of  x  and  y  at  each  point  in  the  field 
between  iterations  is  less  than  the  tolerances  input  as  R(2)  and  R(3) , 
respectively,  or  until  the  maximum  number  of  iterations  allowed  (input 
as  ITER)  is  reached.  In  the  latter  case  the  partially  converged  solu¬ 
tion  is  stored  on  file  10  for  restart.  The  input  parameter  IDISK  can 
cause  the  code  to  read  this  partially  converged  solution  from  file  10 
and  continue  the  iterations.  This  parameter  also  controls  the  dispo¬ 
sition  of  the  final  solution,  which  is  normally  stored  on  file  11  for 
use  in  the  flow  solution,  but  can  be  simply  printed  without  being  stored 
if  desired.  Various  other  input  parameters,  such  as  print  options,  etc., 
are  explained  in  the  detailed  input  instructions  given  below  and  in 
the  source  listing. 

Code  Operation 

Initial  input  and  setup.  The  WESC0R  code  uses  the  values  of  NDIM, 
NDIM1,  NDIM2,  and  NDIM3,  which  are  assigned  by  a  DATA  statement,  to  deter¬ 
mine  if  the  problem  specified  by  the  input  will  fit  in  the  arrays  as 
dimensioned.  The  first  two  of  these  parameters,  NDIM  and  NDIM1,  corre¬ 
spond  to  the  dimensions  of  the  field  arrays,  X,  etc.  The  last  two,  NDIM2 
and  NDIM3,  correspond  to  the  dimensions  of  the  slit  arrays,  XL,  etc.  The 
last  parameter,  NDIM3,  also  corresponds  to  the  dimension  of  the  segment 
arrays,  LB1,  etc.  Thus  NDIM  is  the  maximum  value  of  £  that  can  be  used, 
while  NDIM1  is  the  maximum  value  n  allowed.  Also,  NDIM2  is  the  maximum 
number  of  points  that  can  be  used  on  a  slit  or  slab  side,  and  NDIM3  is 
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Che  maximum  number  of  slits  and  slab  sides  Chat  can  be  used.  The  Input 
thus  must  satisfy  the  following: 

I MAX  <  NDIM 
JMAX  <  NDIM1 

|LB2 (L)  -  LB1 (L)|  +  1  <  NDIM2  L  ■  1,  2,  .  .  . , NBDY 
NBDY  <  NDIM3 

After  the  initial  input  parameters  are  read,  the  code  does  some 
setup  of  various  intermediate  parameters  and  checks  for  comparability 
with  the  array  dimensions.  The  value  of  IDISK  is  then  checked  to  de¬ 
termine  if  the  solution  is  to  be  started  from  the  beginning  or  if  a  par¬ 
tially  converged  solution  is  to  be  continued. 

Boundary  input  and  construction.  If  the  start  is  from  the  beginning, 
the  point  type  array  LSLIT  is  initialized  to  -20000  on  the  outer  rectangle 
formed  by  I  ■  1  &  IMAX  and  J  *  1  &  JMAX,  and  to  0  inside  this  rectangle. 

Next  the  points  on  the  slits  and/or  slab  sides  (if  any)  are  read 
from  either  card  images  or  file  10.  Points  on  slits  are  placed  in  the 
slit  arrays,  XL,  etc.,  while  points  on  slab  sides  are  placed  directly 
in  the  field  arrays  X  and  Y.  The  point  type  array  LSLIT  is  set  to 
-(10000  +  L)  at  points  on  slab  sides,  where  L  Identifies  the  particular 
segment  in  the  order  as  input,  unless  the  side  is  a  part  of  the  outer 
boundary  in  which  case  LSLIT  is  left  at  -20000.  At  the  same  time,  10  is 
added  to  the  segment  type  array  LTYPE  for  slab  sides  on  the  outer  boun¬ 
dary^  resulting  in  replacing  the  input  values  of  -11  and  -12  with  -1  and 
-2,  respectively,  in  conformance  with  the  usage  for  slits. 

The  slit  arrays,  XL,  etc.  (if  any),  are  then  printed  and  subroutine 
BNDRY  is  called  for  the  outer  boundary.  If  the  outer  boundary  is  not 
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input  in  segments  as  slab  sides,  this  boundary  is  either  input  as  a 
succession  of  points  proceeding  from  a  specified  point  completely  around 
the  outer  rectangle  formed  by  I  =  l.IMAX,  and  J  =  1,JMAX,  or  a  circular 
outer  boundary  is  generated  internally  and  placed  on  this  rectangle. 

Both  of  these  procedures  are  performed  by  this  subroutine  by  calling  the 
subroutine  INFBDY,  which  either  reads  a  point  from  a  card  image  or  cal¬ 
culates  a  point  on  the  circle. 


Point  types.  Next  the  point  type  array  LSLIT  is  set  to  the  following 
values  on  and  adjacent  to  slits  (if  any).  Here  L  identifies  the  partic- 


order  as 

input: 

-L  : 

on  slit 

10L  +  1  : 

below  horizontal 

slit 

10L  +  2  : 

above  horizontal 

slit 

not  adjacent 

10L  +  3  : 

left  of  vertical 

slit 

slit  ends 

10L  +  4  : 

right  of  vertical 

slit 

The  point  type  array  LSLIT  is  then  set  to  -10000  for  points  outside 
the  computational  region,  i.e.,  inside  slabs,  by  sweeping  along  each  £- 
and  n-line  and  noting  when  the  computational  region  is  entered  or  left 
across  a  slab  side.  The  complete  point  type  array  then  contains  the 
values  Indicated  in  the  following  diagram: 
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Control  functions  and  Initial  guess.  With  all  of  the  boundary  points 
in  place  and  the  point  type  array  filled,  the  code  then  calls  subroutine 
C0NTRL  to  evaluate  the  control  functions  on  the  entire  boundary  (including 
internal  boundaries).  The  subroutine  GUESSA  is  called  next  to  calculate 
the  control  functions  and  the  initial  guess  for  the  cartesian  coordinates 
in  the  field  by  interpolation  from  the  values  on  the  boundaries.  This 
Interpolation  is  done  at  each  point  in  the  field  by  locating  the  pair  of 
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boundary  segments  facing  the  point  (one  or  both  members  may  be  internal 
boundaries)  and  interpolating  between  these  segments.  For  the  coordinate 
values,  the  distances  separating  the  pair  of  segments  facing  the  point 
in  the  horizontal  and  vertical  directions  are  examined  and  the  interpo¬ 
lation  is  done  between  the  pair  with  the  smaller  separation. 

Iterative  solution.  If  the  solution  is  to  be  restarted  from  a  partially 
converged  result,  then  all  of  the  above  computations  are  skipped  and  the 
partially  converged  solution  is  read  from  file  10  instead.  In  either 
case  the  initial  array  values  are  printed  at  this  point  according  to 
the  input  print  controls. 

Subroutine  TRANS  is  now  called  to  perform  the  iterative  solution. 

This  subroutine  first  reads  the  parameters  associated  with  the  attraction 
of  curvilinear  coordinate  lines  to  other  curvilinear  coordinate  lines  and/ 
or  points.  The  species  of  line  being  controlled,  i.e.,  £  or  n,  is  read 
into  ATYP,  and  whether  the  control  is  to  be  attraction  or  repulsion  is 
determined  by  the  input  parameter  ITYP.  The  number  of  coordinate  lines 
and  points  designated  as  sources  of  attraction  are  read  into  NLN  and  NPT, 
respectively.  Also,  a  common  decay  factor  and  a  common  amplitude  multi¬ 
plication  factor  to  be  used  for  all  attraction  lines  and  points  for  this 
species  can  be  read  into  DEC  and  AMPFAC,  respectively. 

For  each  species  of  control,  subroutine  RHS  is  called  to  read  the 
attraction  line  index,  or  point  indices,  and  the  amplitude  and  decay 
factor  for  each.  This  subroutine  $lso  sums  the  effects  for  all  such 
attraction  lines  and  points  and  adds  this  cumulative  effect  to  the 
control  function  at  each  point  in  the  field  in  accordance  with  Eq.  (5  ) 
of  Part  A. 
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Subroutine  TRANS  then  reads  the  parameters  associated  with  attraction 
of  curvilinear  coordinate  lines  to  specified  lines  and/or  points  in  space 
and  adds  the  cumulative  effect  of  all  such  attraction  lines  and/or  points 
to  the  control  functions  at  each  point  in  the  field.  This  is  done  in  a 
similar  manner  as  described  above.  Subroutine  RHSXY  reads  the  cartesian 
coordinates  of  thepointson  the  specified  attraction  line  and  those  of 
the  isolated  attraction  points  and  calculates  the  normal  to  the  attrac¬ 
tion  line.  These  qualities  are  placed  in  the  arrays  XPT,  YPT,  VEC1,  and 
VEC2.  The  addition  to  the  control  functions  in  this  case  must  be  changed 
as  the  iterative  solution  of  x  and  y  proceeds  since  the  control  functions 
depend  on  x  and  y  for  this  type  of  attraction. 

After  completing  the  calculation  of  the  control  functions,  sub¬ 
routine  TRANS  reads  the  parameters  that  provide  for  a  gradual  implemen¬ 
tation  of  these  equations  during  the  iteration,  and  performs  some  setup 
for  the  iterative  solution. 

The  field  is  then  swept  iteratively  until  convergence  is  achieved 
or  the  maximum  number  of  Iterations  allowed  is  reached.  In  each  itera¬ 
tion,  new  values  for  x  and  y  at  points  having  the  point  type  LSLIT  non- 
negative  are  calculated. 

First,  the  coordinate  derivatives  are  calculated,  and  the  Jacobian 
and  other  such  quantities  and  coefficients  are  evaluated.  Then  the 
locally  optimum  acceleration  parameters  are  calculated  if  such  is  elected. 
The  change  in  these  acceleration  parameters  between  iterations  is  moni¬ 
tored  and  the  values  are  frozen  when  the  magnitude  of  the  change  falls 
below  a  specified  tolerance  at  all  points.  (This  change  between  itera¬ 
tions,  and  the  analogous  changes  in  x  and  y,  are  calculated  by  calling 
subroutine  ERR0R.)  The  acceleration  parameter  is  placed  in  the  field 


65 


iMiMNiteatai 


array  WACC.  The  addition  to  the  control  functions  from  attraction  to 
specified  lines  and/or  points  in  the  physical  region  is  calculated  next, 
and  then  the  new  values  of  x  and  y  for  the  point  are  calculated. 

This  procedure  is  followed  for  all  points  in  the  field,  i.e.,  points 
having  the  point  type  LSLIT  non-negative.  For  points  adjacent  to  slits 
it  is  necessary  to  obtain  the  values  on  the  slit  from  the  slit  arrays, 

XL,  etc.,  and  the  calculations  are  done  in  that  case  by  calling  sub¬ 
routine  SLIT. 

After  each  sweep  of  the  field  the  maximum  changes  in  x  and  y  from 
the  previous  sweep  are  compared  with  the  input  tolerances.  If  the  max¬ 
imum  number  of  iterations  allowed  by  the  input  is  reached  before  conver¬ 
gence,  then  the  partially  converged  solution  is  written  on  file  10  for 
potential  restart.  If  convergence  is  obtained,  the  solution  is  written 
on  file  11. 

LINES  (Boundary  Segments) 

The  small  front-end  code  LINES  generates  a  distribution  of  a  speci¬ 
fied  number  of  points  on  a  curve  between  two  specified  points.  The  curve 
may  be  specified  to  be  a  straight  line,  a  circular  or  elliptic  arc,  a 
quadratic  with  zero  slope  at  either  end  point,  or  a  cubic  with  the  slope 
specified  at  both  ends.  In  any  case  the  point  distribution  on  the  curve 
may  be  uniform  or  exponentially  concentrated  toward  either  end.  The 
input  consists  of  the  number  of  curves  to  be  generated  and,  for  each 
curve,  the  number  of  points  on  the  curve,  the  type  of  curve,  the  end 
points,  and  the  particular  quantities  to  be  specified  in  connection  with 
each  curve.  Detailed  instructions  for  input  are  given  below. 
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The  cartesian  coordinates  of  the  points  generated  on  each  curve 
are  output  In  succession  on  file  10  by  a  separate  unformatted  write 
statement  for  each  point  (WRITE(IO)  X,Y).  Since  more  than  one  curve 
can  be  generated  In  one  run,  this  code  can  be  used  to  build  an  entire 
boundary  composed  of  segments  of  different  types.  The  generation  of  the 
curves  and  the  exponential  concentration  of  points  thereon  are  explained 
in  the  following  section. 

Generation  of  Curves 

/ 

Straight  line.  Here  we  have  simply 

y  *  a  +  bx 

so  that  with  the  end  points  (x^,  y^)  and  (xj,  y2)  specified  we  have 


1  x2 

V 

m 

’i 

1  *2, 

lbi 

r* 

so  that 


ylX2  "  y2Xl 


Circular  arc.  For  a  circular  arc  of  radius  r  centered  at  (x^,  y^) 
with  6  measured  counter-clockwise  from  the  positive  x-axis,  we  have 
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x  *  Xq  +  r  cob  6 
y  -  yQ  +  r  sin  0 

The  end  points  are  defined  by  inputting  the  radius  r  and  center  of 
the  arc  (Xq,  yg).  together  with  the  angles  0^  and  @2  of  the  end  points. 


Elliptic  arc.  In  this  case  we  have,  for  an  ellipse  with  semi-major 
axis,  a,  and  semi-minor  axis,  b,  centered  at  xQ,  yQ,  the  equation 


(x  -  xQ)2  (y  -  yQ)2 
_2 - + - ^ - 


which  can  be  written  in  terms  of  the  angle  0,  measured  counter-clockwise 
from  the  positive  x-axis,  and  the  angular-dependent  radius  r(0)  as 

x  ■  Xq  +  r(0)  cos  0 

y  -  yQ  +  r(0)  sin  0 


Then 


r(0) 


cos2  0  .  sin2  0  ^ 

“ ^“  +  “b2— 
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i  end  points  are  specified  by  inputting  the  axes  a  and  b,  the  center 
l<  Yq)»  and  the  angles  of  the  end  points. 


Quadratic  with  zero  slope  at  end  point.  Here  we  have 

y  ■  a  +  bx  +  cx2 
y "  =  b  +  2cx 

Then  with  the  end  points  (x^  y1)  and  (x2>  y2>  specified  together 
with  the  specification  of  zero  slope  at  end  point  i  (i  *  1  or  2)  we  have 


which  is  solved  for  the  coefficients  a,  b,  c. 
Cubic.  The  cubic  equation  is 

y  -  a  +  bx  +  cx2  +  dx3 

y'  -  b  +  2cx  +  3dx2 


or 
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which  is  solved  for  the  coefficients  a,  b,  c,  d. 

Exponential  Concentration  of  Points 

The  exponential  distribution  of  points  on  the  curve  of  any  type  is 

done  by  taking 


for  concentration  near  the  first  end  point  and 


for  concentration  near  the  second  end  point.  Here  the  strength  of  the 
concentration  is  controlled  by  the  specified  decay  factor  a,  and  N  is 
the  number  of  points  on  the  curve. 
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CSPL0T  (Plot) 

The  plot  code  CSPL0T  plots  the  coordinate  system  generated  by  the 
code  WESC0R,  having  read  the  coordinate  system  from  file  11  as  output 
by  WESC0R.  The  Input  consists  of  the  number  of  coordinate  lines  to  be 
plotted,  a  designation  for  skipping  lines,  the  extent  of  the  field  to  be 
plotted,  and  a  factor  for  using  different  seating  in  the  horizontal 
and  vertical  directions.  This  input  is  detailed  in  the  following  pages. 

The  plot  Is  formed  by  simply  connecting  the  points  on  a  line  of 
constant  curvilinear  coordinates  in  the  physical  region,  i.e.,  by  con¬ 
structing  straight  lines  between  each  successive  pair  of  points,  X(I,J) 
and  Y(I,J),  as  one  index  is  held  fixed. 
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WESC0R  INPUT  INSTRUCTIONS 


saw 

m-t 

w*c 

31<K 

«20=C 

°30=C 

caO-C 

9tHK 
”0=C 
®do=.; 
1»0=C 
lOOO=C 
1J1G=C 
1020=C 
10 30=C 
:040=c 
1«0=C 
1060-C 
1070=C 
1080=C 
:09o=c 
uoo=c 
:iiO=C 
:i20=C 
1'30=C 
;i40=c 

uso=c 

ilftO-C 
U70--C 
:i80=C 
i:90=c 
i200=C 
1210- C 
1220=C 
1230=C 
1240=C 
1330=C 
1260=0 
1270=C 

i:-8o=c 

1290=C 

1300=C 

1310-C 

1320=C 

1330=C 

1340=0 

1350=C 

1340=0 

1370=0 

1380=0 

1390=0 

1.00=C 

1410-C 

1420=C 

1430=0 

1140=0 

1450=0 

1160=C 

1470=0 

li80=C 

M90=C 

1500=0 

1510=0 

i520=C 

1530=0 

1540=0 

1550=C 

1500=0 

1570=0 

1580=0 

1590=0 

lo00=C 

1610=C 

lo20=C 

1>30=C 

lo40=C 

U50=C 


"REN  -  NGN-ZERO  USES  RE-EnTWWC  BOUNDARY  08  L£FT  4  RICrlT  SUES 
OF  TRANSFORMED  RECI3M.  WITH  OUTER  BOUNDARY  ON  TOP 

and  inner  boundary  on  bottom. 

inner  boundary  is  road  as  rollons  before  reading  outer  ioumart: 

=1  INNER  BOUNDART  READ  FRON  CARDS. 

X.Y  -  FORMAT!  2f»0.0>  .  ONE  POINT  PER  CARD. 

=2  inner  boundary  read  from  file  ig. 

X.Y  -  UNFORMATTED  .  ONE  INAGE  PER  CARD. 

(NOTE!  SLITS  AND/OR  SlAIS  NAT  AtSO  K  PRESENT.) 

F0RNATmI5) 


m  CARDS!  NBBY)  L31.LB2.LB3.LIYPE 

t 

iBl.LB.2  -  FIRST  AnD  LAST  INDICES  OF  kAB  SIDE  OR  SlIT  ENDS. 

(LD2  NAT  DC  LESS  TnAM  LSI  FOR  SLAB  SIDE.  INPUT  IS  FRON  LB1  TO  LD2. ) 


I 
t 
I 
I 
I 
t 
« 
t 
> 

< 
t 
t 

III  CARD 

« 
t 
I 
t 
t 
I 
t 
I 
I 
I 
t 
* 

I 

« 

* 

I 
t 
« 

1 
I 
I 
« 
t 

I . 

I 

U  IF  BODIES  AND/OR  OUTER  BOUNDARY  ARE  READ  FRON  CARDS’  SUCH  CARDS 

II  FOLLOW  NEXT. 

II 

It  SilTS  AnD/OR  Si. AD  SIDES  ARE  READ  FIRST.  THEN  OUTER  BOUNDARY  IS  READ. 

I  .  THESE  RULES  APPLY  FOR  READING  FRON  FILE  10  AS  DELL  AS  FRON  CARDS. ) 

II 

I . 

I 

II  If  (.0  COORDINATE  ATTRACTION  IS  TO  BE  USED.  EOLLON  THESE  CARDS 

II  WITH  FIVE  Blank  CARDS.  If  ATTRACTION  IS  TO  BE  USED.  USE  THE  FOLLOWING 

it  input  rather  than  ThE  blank  cards: 

It 

u  input  for  coordinate  systen  control  :  use  four  sets,  one  for 

II  XI-LINE  ATTRACTION  TO  COORDINATE  UNTS/fOINTS.  ONE  FOR  ETA-LInE  ATTRACTION 
II  TO  COORDINATE  UNE3/P0INTS.  ONE  FOR’  XI -LINE  ATTRACTION  TO  SPACE  LINES/POINTS. 
II  AnB  ONE  FOR  ETA-UNE  ATTRACTION  TO  SPACE  LINES/ROINTS. 

U  ANY  SET  NOT  RANTED  IS  REPLACED  BY  ONE  BLANK  CARD. 

II 

tHKmTmHIttmmilimtlHMIMtHHMmiMHmHBMIlNI 

I* 

II  THE  FOLLOWING.  NARKED  WITH  I-  IS  FOR  ATTRACTION  TO  COORDINATE  LINES/POINTS: 

II 

1411  CARD  :  ATYFi!TYPiir..N.Nf'TiBfC.ANRfAC  -  FORNATl  A8.I2.2I5.2E1C.0) 

II 

II  A  TIP  -  TYPE  OF  ATTRACTION.  Ui  FOR  XKINE  ATTRACTION. 

II  ETA  FOR  ETAHINE  ATTRACTION.)  LEFT  JUSTIFIED. 


LBS  -  INDEX  OF  LINE  On  WHICH  Si. AD  SIDE  OR  S.IT  IS  LOCATED. 

LTTPE  -  SLAB  SIDE  OR  Si.IT  TYPE  (1  FOR  HORIZONTAL.  2  FOR  VERTICAL.) 

(negative  indicates  slab  side,  rather  than  slit.) 

i SUBTRACT  10  FOR  OUTER  BOUNDART  SEGNENT.  ) 

(I.F..  -11  IS  HORIZONTAL  OUTER  BOUNDARY  SEGNENT.) 

(  -12  IS  VERTICAL  OUTER  BOUNDARY  SEGNENT.  I 

R<  1  "Ri2'fR(3»YlNf!N>Alff  1N.X0INF .  YOIWF.  INFXI  •  InFETA 
-  FORMAT*  7F10.0.2I5) 

R(l)  -  SOR  ACCELERATION  RARANETER. 

(ZERO  VAluE  CAUSES  VARIABLE  ACCELERATION  PARAMETER) 

(FIElD  TO  BE  CALCULATED  INTERNAlLY.  ) 

Rt 2 )  -  ALLOMABi.E  X  ITERATION  ERROR. 

R(3)  -  ALLOWABLE  Y  ITERATION  ERROR'. 

YINFIN  -  RADIUS  OF  CIRCULAR  OUTER  BOUNDARY. 

AINFIN  -  ANGLE  OF  FIRST  ROINT  ON  CIRCULAR  OUTER  10UNDARY  (DEGREES). 

( COUNTER-CLOCK  FRON  POSITIVE  X-AXIS.) 

XOINF.YOINF  -  CENTER  OF  CIRCULAR  OUTER  BOUNDARY. 

NINF  •  NUMBER  OF  UNIDUE  POINTS  On  CIRCULAR  OUTER  BOUNDARY. 

ihfxi. infeta  -  indices  of  first  point  on  circular  outer  boundary. 

(note  :  LAST  7  OF  TnESE  PARAMETERS  ARE  IRRELEVANT  IF  OUTER  BOUNDARY  IS  READ.) 
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TIC  LAST  COORDINATE  STSTEN  CONTROL  CARD  IS  THE  FOLLOWING  CARD  ! 

CARD  !  IFAC.IRITiEFAC  -  FORMAN  215.F10.0) 

(CAN  BE  USED  TO  AID  CONVERGENCE  BY  CONVERGING  FIELD  ) 

(WITH  LESS  ATTRACTION  FIRST  AND  USING  THIS  RESULT  ) 

(AS  THE  INITIAL  GUESS  FOR  STRONGER  ATTRACTION.  ) 
(BLANK  CARD  MUST  BE  INPUT  IF  THIS  FEATURE  IS  NOT  USED. ) 
(STANDARD  IS  TO  NOT  USE  THIS  FEATURE  >  BUT  ITS  USE  NAY) 

(BE  NECESSARY  KITH  STRONG  ATTRACTION.  ) 

IFAC  -  NUMBER  OF  STEFS  IN  ADDITION  OF  INHOMOGENEOUS  TERM. 
DOUBLES  INHOMOGENEOUS  TERN  AT  EACH  STEP. 

(ZERO  CONVERGES  WITH  FULL  ATTRACTION 
<1.0  CONVERGES  MITH  NO  ATTRACTION  FIRST •  TtCN 
(HITH  FULL  ATTRACTION.  2.0  CONVERGES  WITH  NO  > 
(ATTRACTION  FIRST.  THEN  HITH  HALF.  THEN  WITH  FULL.) 

( INCREASE  NUMBER  OF  STEFS  IF  DIVERGENCE  OCCURS.  > 

IRIT  -  NON-ZERO  VAlUE  CAUSES  INHOMOGENEOUS  TERM  TO  BE  PRINTED. 

EFAC  -  MULTIPLE  OF  CONVERGENCE  CRITERION  TO  BE  USED  FOR 
INTERMEDIATE  CONVERGENCE  BETWEEN  ADDITIONS  OF 
INHOMOGENEOUS  TERM.  ( TYPICALLY  10.0  ) 


NASS  STORAGE  FILES  ! 
RESTART  FILE  -  FILE  10  I 


(10)  RXI.RETA 

(10)  X.Y.LSLIT.LABEL.IMAXiJMAX 
(10)  NBDY.NUHB.LB1.LB2.LB3.LTYPE.LPT.XL.XU.YL. YU. 
NDIH.NDIH1.NBIN2.NDIH3.UACC 


COORDINATE  SYSTEM  STORAGE  FILE  -  FILE  11 


(11)  LABEL.IMAX.JHAX 

(11)  (USLIT(I.J)>]-1>]MAX)>J=1..1MAX) 

(ID  (( X(I.J).  1=1  .IMAX).J=1.  UMAX) 

(11)  <(Y(I«J)rI=li IHAX )r J=1  •  JHAX ) 

( 11 )  NBDYiNUMBiLB1iLB2iLB3iLTYPE.LPT.XL.XUiYL.YU. 
NDIH.NDIM1.ND1H2.ND1M3 


LINES  INPUT  INSTRUCTIONS 


:40=C« 

;io=c 

120=C* 

130=C 


LIMES  wmuamuwmnmmmiui 


'S0=C 

i«0=C  BOUNDARY  SEGMENT  CODE  FOR  INPUT  TO  MESCOfi 
170=C 

180=C  MISSISSIPPI  STATE  UNIVERSITY  .  19S2 

:*o=c 

200=C  U.S.ARAY  ENGINEER  UATERUAYS  EXPERIMENT  STATION 

:iO=C  VICKSBURG,  MISSISSIPPI 

224=C 

'Mmtimmtmnmuiimmmmimmimmummmmmtt 

140=5 

bo=COT*  POINTS  ON  BOUNDARY  SEGMENTS  OSS 

2ao=c 


290=ctm  input  : 

3g^C 

3NKIICARB  !  KLINES  -  FORMAT!  15) 

J20=C 

130*C  KLINES  -  TOTAL  NUHIER  OF  LINES. 

340=5 

TSO^CtSCARDSf KLINES )  !  N,ITYP,Dl,D2,B3,B4,B5,Il6,DF.  - 
360=0 

N  -  NUNBER  OF  POINTS  ON  LINE. 


FORMAT! 2I5.7F10.0) 


37<K 

380=C 

390=C 

400=0 

*10=C 

♦20=C 

•30=C 

♦40=C 

*50=C 

»40=C 

»70=C 

*80=C 

*90=C 

500=C 

S10=C 

520=C 

530=C 

540=C 

550=C 

540=C 

570=C 

580=C 

390=C 

sOO=C 

410=C 

o20=C 

430=C 

o40=C 

o50=C 

o40=C 

470=0 

«80=C 

«?0=C 

700=0 

710=0 

720=0 

730=5 

740=0 

TO=C 

740=0 

770=C 

780=C 

790=0 


ITYP  -  TYPE  OF  LINE 
0  i  STRAIGHT. 

2  :  CIRCULAR  ARC. 

2  !  ELLIPTIC  ARC. 

3  !  CUBIC. 

4  :  QUADRATIC  UITH  ZERO  SLOPE  AT  FIRST  POINT. 

5  !  QUADRATIC  UITH  ZERO  SLOPE  AT  SECOND  POINT. 

Dl-Dt  AS  FOLLOWS  -  ( ITEMS  NOT  CITED  ARE  IRRELEVANT) 

ITYP=0  !  D1  -  X  OF  FIRST  POINT. 

D2  -  T  OF  FIRST  POINT. 

B3  -  X  OF  SECOND  POINT. 

D4  -  Y  OF  SECOND  POINT. 

ITYP=1  :  III  -  ANGLE  OF  FIRST  POINT  (DEGREES. 

D2  -  ANCLE  OF  SECOND  POINT  (DEGREES. 

D3  -  X  OF  CIRCLE  CENTER. 

D4  -  Y  OF  CIRCLE  CENTER. 

D5  -  CIRCLE  RADIUS. 

ITYP=2  )  Dl  -  ANGLE  OF  FIRST  POINT.  ( DECREES. 

D2  -  ANGLE  OF  S£C»  POINT.  (DEGREES, 

D3  -  X  OF  EU.IPSE  CENTER. 

D4  -  Y  OF  ELLIPSE  CENTER. 

DS  -  X-AXIS  LENGTH  OF  ELLIPSE. 

D6  -  Y-AXIS  LENGTH  OF  aLIPSE. 

ITYP=3  5  Dl-Di  SANE  AS  ITYP=0 

D5  -  SLOPE  AT  FIRST  POINT.  (DECREES, 

D4  -  SLOPE  AT  SECOND  POINT.  (DEGREES, 


COUNTER -CLOCK  FROM  POSITIVE  X-AXIS) 
COUNTER-CLOCK  FROM  POSITIVE  X-AXIS) 


COUNTER-CLOCK  FROM  POSITIVE  X-AXIS) 
COUNTER-CLOCK  FROM  POSITIVE  X-AXIS) 


COUNTER-CLOCK  FROM  POSITIVE  X-AXIS) 
COUNTER-CLOCK  FROM  POSITIVE  X-AXIS) 


ITYP*4  5  D1-D4  SANE  AS  ITYP=« 

ITYP=5  J  D1-D4  SAME  AS  ITYf=4 

K  -  EXPONENTIAL  CONCENTRATION  FACTOR. 

4.0  FOR  EQUAL  SPACING  ON  LINE. 

NEGATIVE  FOR  CONCENTRATION  NEAR  FIRST  POINT, 
POSITIVE  FOR  CONCENTRATION  NEAR  SECOND  POINT. 


D14=C 

820=0  MASS  STORAGE  FILE 
88  OUTPUT  -  FILE  10  5 

tm.p 

wrv 


URITE(IO)  X(I>,Y(I> 

X  4  Y  POINTS  OF  EACH  LINE, 


INCLUDING  ENDS. 
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CSPL0T  INPUT  INSTRUCTIONS 


SAMPLE  RUNSTREAMS 


LINES  Sample  Runstream  #1 


120= JOB.  _ 

:»«4CaiRE.^BINARt»PW*TH(W>SmiNESJ.IIl4BB*I.Iir=TRiUfl. 

140=LJR,ttf=BlNARY .  _ 

150=BISPOSE,BN=FT10,SBN=FILE>  ID=MM»>  KfSTrBF=TR»TEXT=t 


loO=  'CATALOG,  FILE,TH0NPS0NUN£S0,Ii 
170=0£LETE*BN=BINARY . 

180=EXIT. 

lPO=DaETE  >  D#=8IHARY . 

2oo=« 


27O=*£0R 

28O=t£0F 


»RP=999. 


210= 

220= 

5 

33 

0  0.0 

0.0 

24.39 

0.0 

230= 

9 

0  0.0 

-0.3 

6.1 

-0.3 

240= 

25 

0  A.l 

-0.3 

24.39 

-0.91 

250= 

25 

0  0.0 

-0.3 

0.0 

0.0 

260= 

25 

0  24.39 

-0.91 

24.39 

0.0 

WESC0R  Sample  Runstream  //I 


120=J(»,T=60.  _ 

l30=ACflUIRE,HH:T10,PIIAgTH0ltPS0WLlllES0ilD=WWIrI|E=Tll,U0. 
140=AC#1IRE ,  BK=BIMARY ,  PBM=THQHPSOrtCCIRBB » IMHil,  DF=TR , UQ . 

150=UR,DN=81I1ARY.  _ 

loO=DISfOSE,DH=fTll.SDM=FIL£,ID=«W —  »DC=ST>HF=TR»TEXT=t 
170=  'CATAL0G,FILE,TH0HPS0MC0RB0,ID=flBBBiRP=999.' . 
180= DELETE ,  WP=FT 1 0 . 

190=  DELETE (ORDINARY. 

20o=exit. 

210=B£LETE,D#=FT10. 

220=B£LETE,D*=BINARY. 

23O=tE0R 

24O=JBKS0N  FLUME 
250*33  X  25  COORDINATE  SYSTEM 


2»0= 

270= 

280= 

290= 

300= 

310= 
320=1.8 
j30= 

340= 

350* 

3e0= 

370= 

380=»E0R 

390=*£0f 


25 

33 

9 

33 

25 

25 


100 

-11 

-11 

-11 

-12 

-12 


-1 


1  1 


0 


0.00001  0.00001 


CSPL0T  Sample  Runstream 


i,DF=TR,U0. 
DF=TR,UQ. 


120*001 

l30®AOSJIRE»DN=fTlO,PBN=THONPSONCORBO,ID: 
140=AMUIRE,BM*BIIW(Y,PDN=THOMPSONCSPLOTBri; 

150*LIR,LIB=METALIB,BN=BIIIARY.  _ 

1AO=DISPOSE,BM*ET01,SM<=FILE,1^1— l,DC=ST,DF=SB,T£XT=t 
170=  'CATAt.K,FILE,THOMPSOMPLOT,ID«MMI,RP=9?9.'. 
180= DELETE, DM*FT10. 

190*D£LETE,DM*»INARY. 

200=EXIT. 

2lO=DELETE,DN*fT10. 

220* DELETE, BM»IIMARY. 

230=K» 

240*  0  0  1  1 

250*  0  0  0  0 

200=2.0 
27C 
210*1 
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Bi 


LINES  Sample  Runstream  j 


;20=J(*.  _____ 

1 30=AC0UI  kE » iH=B  I  HAST  >  PDM= THOMPSOHL  1HES&  r  I  I>r  =TR  >  UQ 

140^K'M>tINA*Y. 

t50=»lSW$E ,  Bm=FT  10>SIlM=f  ILEfI  DHM  >HC=ST>DF=Tt<.TEXT=t 
140=  'UTALOG,FlLE.Th6rtPmiNES01.IO=*B*»,Rf=999,' 
170=IEL£TE»M=BI)M*:Y. 


160=£XIT. 

190=KLETE,M=KNART. 

200=lfOR 


:io= 

15 

v>0= 

25 

0  15.0 

0.0 

:so= 

9 

0  5.91 

0.0 

340= 

3 

0  4,1 

0.48 

250= 

9 

0  3.9 

0.48 

:«o= 

7 

0  2.29 

0.0 

’70= 

*n 

0  0.0 

0.0 

280= 

21 

0  0.0 

1.43 

2°0= 

29 

0  5.91 

1.43 

300= 

7 

0  15.0 

1.43 

310= 

5 

3  15.0 

1.25 

320= 

5 

0  15.25 

330= 

0  15.8 

nwVj 

340= 

0  15.43 

■If/I 

350= 

3  15.25 

0.58 

jo0= 

ro=«E« 

380= It Of 


5.91 

0.0 

0.0 

4.1 

0.48 

0.0 

3.9 

0.48 

2.29 

0.0 

0.0 

0,0 

0.0 

0,0 

0.0 

l.o3 

0.0 

5.91 

1.43 

0.0 

15.0 

i  .43 

0.0 

15.0 

1.25 

15.25 

0.94 

-57.0 

15.8 

0.43 

0.0 

15.43 

0.23 

0.0 

1 '  4 

0.58 

0.0 

1-  ■ 

HH 

133.0 

Ij  ■ 

■flMl 

0.0 

WESC0R  Sample  Runstream  //2 


120=J08»T=60. 

l34=«CflUIRE>ltt=FT10.PM=TH0llf’S0«lIftES01>lIF 
140=MCQUlR£>IM=BINtftYrPM=TM)HPS(MCOKI)8>II|: 

150=114  iDM*HNA8Y  > 
l&0=IISro&EfM=ETU>SM=fILE>ni=««HI>DC=ST>tF=TRiT£XT=t 
1/0=  'CAIAL0t>FIUMhOYS0C0ftIl01>I&  " 

180=Kl£TE>M=FT10. 
l90=KLETE.W*=HMtfT. 

200=£XJT. 

210=ICLETE»W=fT10. 

220=Ki.ETE  iM=81M8Y . 

230=«* 


iDf=TE.U6. 

>Df=TK.U0. 


MP=99?.\ 


744=  MKTCM  TEST  4  -  KITH  HEIR 
250=  57  X  23  COQMINMTE  SYSTEM 

240= 

57 

23 

15  100  2 

270= 

4* 

25 

1  -11 

280= 

25 

17 

1  -11 

290= 

i7 

17, 

1  -11 

300= 

15 

7 

1  -11 

310= 

7 

1 

1  -11 

320= 

i 

23 

1  -12 

330= 

i 

7; 

23  -11 

340= 

21 

49 

23  -11 

350= 

23 

17 

49  -2 

j«0= 

49 

53 

1/  -1 

370= 

33 

57 

17  -1 

380= 

17 

9 

57  -12 

390= 

57 

53 

9  -j 

400  = 

53 

4? 

9  -1 

410= 

9 

1 

49  -2 

420=1 

.8 

0.00001  0.00001 

-1 


1  1 


430= 
440= 

450= 

4«0= 

470= 

4*0=*t« 

49O=lE0f 


(Reduced  Horizontal 


Scale) 


Penstock 


Weir 


In  accordance  with  letter  from  DAEN-RDC ,  DAEN-ASI  dated 
22  July  1977,  Subject:  Facsimile  Catalog  Cards  for 
Laboratory  Technical  Publications,  a  facsimile  catalog 
card  in  Library  of  Congress  MARC  format  is  reproduced 
below. 


Thompson,  Joe  F. 

A  boundary-fitted  coordinate  code  for  general  two- 
dimensional  regions  with  obstacles  and  boundary 
intrusions  /  by  Joe  F.  Thompson  (Department  of  Aerospace 
Engineering,  Mississippi  State  University).  --  Vicksburg, 
Miss.  :  U.S.  Army  Engineer  Waterways  Experiment  Station  ; 
Springfield,  Va.  ;  available  from  NTIS,  1983. 

80  p.  :  ill.  ;  27  cm.  --  (Technical  report  ;  E-83-8) 

Cover  title. 

"March  1983." 

Final  report. 

"Prepared  for  Office,  Chief  of  Engineers,  II  S.  Army 
under  EWQOS  Task  IIIA.4." 

"Monitored  by  Hydraulics  Laboratory,  U.S.  Army  Engineer 
Waterways  Experiment  Station." 

At  head  of  title:  Environmental  5  Water  Quality 
Operational  Studies. 

Bibliography:  p.  48. 


Thompson,  Joe  F. 

A  boundary-fitted  coordinate  code  for  general  :  ...  1983 

(Card  2) 

1.  Boundary  value  problems.  2.  Computer  programs. 

3.  Coordinates.  4.  Difference  equations.  5.  Differential 
equations.  6.  Finite  difference  equations.  7.  Numerical 
analysis.  I.  Mississippi  State  University.  II.  United 
States.  Army.  Corps  of  Engineers.  Office  of  the  Chief 
of  Engineers.  III.  Environmental  5  Water  Quality 
Operational  Studies.  IV.  U.S.  Army  Engineer  Waterways 
Experiment  Station.  Hydraulics  Laboratory.  V.  Title 
VI.  Series:  Technical  report  (U.S.  Army  Engineer  Waterways 
Experiment  Station)  ;  E-83-8. 

TA7.W34  no. E-83-8 


